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17 Abstract. 

is Further development of our self-consistent model of interacting ring current 

19 (RC) ions and electromagnetic ion cyclotron (EMIC) waves is presented. This model 

20 incorporates large scale magnetosphere-ionosphere coupling and treats self-consistently 

21 not only EMIC waves and RC ions, but also the magnetospheric electric held, RC, 

22 and plasmasphere. Initial simulations indicate that the region beyond geostationary 

23 orbit should be included in the simulation of the magnetosphere-ionosphere coupling. 

24 Additionally, a self-consistent description, based on first principles, of the ionospheric 

25 conductance is required. These initial simulations further show that in order to model 

26 the EMIC wave distribution and wave spectral properties accurately, the plasmasphere 

27 should also be simulated self-consistently, since its fine structure requires as much 

28 care as that of the RC. Finally, an effect of the finite time needed to reestablish a 

29 new potential pattern throughout the ionosphere and to communicate between the 

30 ionosphere and the equatorial magnetosphere cannot be ignored. 
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si 1. Introduction 

32 Electromagnetic ion cyclotron (EMIC) waves are a common and important feature 

33 of the Earth’s magnetosphere. The source of free energy for wave excitation is provided 

34 by the temperature anisotropy of ring current (RC) ions, which naturally develops during 

35 inward convection from the plasmasheet. The EMIC waves have frequencies below the 

36 proton gyro-frequency, and they are excited mainly in the vicinity of the magnetic 

37 equator with a quasi field-aligned wave normal angle [ Cornwall , 1965; Kennel and 

38 Petschek, 1966]. These waves were observed in the inner [ LaBelle et al, 1988; Erlandson 

39 and Ukhorskiy, 2001] and outer [Anderson et al, 1992a, 1992b] magnetosphere, at 

40 geostationary orbit [Young et al, 1981; Mauk, 1982], at high latitudes [Erlandson et al, 

41 1990], and at ionospheric altitudes [Iyemori and Hayashi, 1989; Braysy et al, 1998]. 

42 Feedback from EMIC waves causes nonadiabatic pitch-angle scattering of the RC 

43 ions (mainly protons) and their loss to the atmosphere, which leads to the decay of RC 

44 [Cornwall et al, 1970]. This is especially important during the main phase of storms, 

45 when RC decay is possible with a time scale of around an hour or less [Gonzalez et al, 

46 1989]. During the main phase of major storms RC 0 + may dominate [Hamilton et al, 

47 1988; Daglis, 1997]. These ions cause damping of the T/e + -mode EMIC waves, which 

48 may be very important for RC evolution during the main phase of the greatest storms 

49 [Thorne and Horne, 1994; 1997]. Obliquely propagating EMIC waves interact well with 
so thermal plasmaspheric electrons due to Landau resonance [Thorne and Horne, 1992; 

5 i Khazanov et al, 2007b], Subsequent transport of the dissipating wave energy into the 
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52 ionosphere causes an ionospheric temperature enhancement [ Gurgiolo et al, 2005]. This 

53 wave dissipation is a mechanism proposed to explain stable auroral red arc emissions 

54 present during the recovery phase of storms [Cornwall et al., 1971; Kozyra et al., 1997]. 

55 Measurements taken aboard the Prognoz satellites revealed a so-called “hot zone” 

56 near the plasmapause, where a temperature of plasmaspheric ions can reach tens of 

57 thousands of degrees [ Bezrukikh and Gringauz, 1976; Gringauz, 1983; 1985]. Nonlinear 

58 induced scattering of EMIC waves by thermal protons [ Galeev , 1975] was used in the 

59 RC-plasmasphere interaction model by Gorbachev et al. [1992] in order to account for 
eo these observations. An extended analysis of thermal/suprathermal ion heating by EMIC 
ei waves in the outer magnetosphere was presented by Anderson arid Fuselier [1994], 

62 Fuselier and Anderson [1996] and Horne and Thorne [1997]. Relativistic electrons 

63 (> 1 MeV) in the outer radiation belt can also strongly interact with EMIC waves 

64 [ Thorne and Kennel, 1971; Lyons and Thorne, 1972], Data from balloon-borne X-ray 

65 instruments provides indirect but strong evidence that EMIC waves cause precipitation 
ee of outer-zone relativistic electrons [ Foat et al., 1998; Lorentzen et al, 2000]. These 

67 observations stimulated theoretical and statistical studies, which demonstrated that 
es EMIC wave-induced pitch-angle diffusion of MeV electrons can operate in the strong 

69 diffusion limit with a time scale of several hours to a day [Summers and Thorne, 2003; 

70 Albert, 2003; Meredith et al, 2003]. This scattering mechanism is now considered to 

71 be one of the most important means for relativistic electron loss during the initial and 

72 main phases of storm. All of the above clearly demonstrates that EMIC waves strongly 

73 interact with electrons and ions of energies ranging from ~ 1 eV to ~ 10 MeV, and that 
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74 these waves strongly affect the dynamics of resonant RC ions, thermal electrons and 

75 ions, and the outer radiation belt relativistic electrons. The effect of these interactions is 

76 nonadiabatic particle heating and/or pitch-angle scattering, and loss to the atmosphere. 

77 The rate of ion and electron scattering/heating in the Earth’s magnetosphere is not 

78 only controlled by the wave intensity-spatial-temporal distribution but also strongly 

79 depends on the spectral distribution of the wave power. Unfortunately, there are still 
so very few satellite-based studies of EMIC waves, especially during the main phase of 

81 magnetic storms, and currently available observational information regarding EMIC 

82 wave power spectral density (mainly from the AMPTE/CCE and CRRES satellites) 

83 is poor [Engebretson et al, 2008 ]. Ideally, a combination of theoretical models and 

84 available-reliable data should be utilized to obtain the power spectral density of EMIC 

85 waves on a global magnetospheric scale throughout the different storm phases. To 

86 the best of our knowledge, there is only one model that is able to self-consistently 

87 simulate a spatial, temporal and spectral distribution of EMIC waves on a global 

88 magnetospheric scale during the different storm phases [ Khazanov et al., 2006 ]. This 

89 model is based on first principles, and explicitly includes the wave generation/damping, 

90 propagation, refraction, reflection and tunneling in a multi-ion magnetospheric plasma. 

91 The iJe + -mode EMIC wave simulations based on this model have showed that the 

92 equatorial wave normal angles can be distributed in the source region, i. e. in the region 

93 of small wave normal angles, and also in the entire wave region, including those near 

94 90 °. The occurrences of the oblique and field-aligned wave normal angle distributions 

95 appear to be nearly equal with a slight dominance of oblique events [Khazanov and 
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ge Gamayunov, 2007]. This theoretical prediction is supported by a large data set of the 

97 observed wave ellipticity [Anderson et al, 1992b; Fraser and Nguyen, 2001; Meredith et 

98 al., 2003]. The observation of a significant number of linearly polarized events near the 

99 equator suggests that waves are often highly oblique there. Using the more reliable wave 

100 step polarization technique, Anderson et al. [1996] and Denton et al. [1996] analyzed 

101 data from the AMPTE/CCE spacecraft, presented the first analysis of near linearly 

102 polarized waves for which the polarization properties were determined. They found a 

103 significant number of wave intervals with a wave normal angle 6 > 70°, the highest 6 

104 ever reported. Compared to field-aligned waves, such highly oblique wave normal angle 

105 distributions can dramatically change the effectiveness (by an order of magnitude or 
loe more) of both the wave-induced RC proton precipitation [Khazanov et al., 2007b] and 
107 relativistic electron scattering [Glauert and Horne, 2005; Khazanov and Gamayunov, 
loe 2007]. Strong sensitivity of the scattering rates to the wave spectral characteristics, 

109 and the wide distribution of EMIC wave normal angles observed in the magnetosphere, 
no suggests that in order to employ EMIC waves for heating and/or scattering of the 

in magnetospheric particles in a model, the wave spectral distribution will require special 

112 care, and should be properly established. 

113 The resulting EMIC wave power spectral density depends on the RC and cold 

114 plasma characteristics. On the other hand, the convective patterns of both RC ions 

115 and the cold plasmaspheric plasma are controlled by the magnetospheric electric field, 
no determining the conditions for the interaction of RC and EMIC waves. Therefore, this 
m electric field is one of the most crucial elements necessary to properly determine the 
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ns wave power spectral density. The region 2 field-aligned currents (FACs) couple the 

119 magnetosphere and ionosphere. This large scale coupling determines and maintains a 

120 self-consistent dynamic of the electric field and RC [Vasyliunas, 1970; Jaggi and Wolf, 

121 1973; Garner et al, 2004; Fok et al, 2001; Khazanov et al, 2003b; Liemohn et al, 

122 2004], A self-consistent simulation of the magnetosphere-ionosphere system should 

123 provide, at least in principle, the most accurate theoretical electric field. The EMIC 

124 waves resulting in the magnetosphere are not only a passive element in the coupled 

125 RC-ionosphere system but also may influence the electrodynamics of coupling. During 

126 storm times, the wave-induced RC proton precipitation not only changes the FAC 

127 distribution, but can potentially modify the conductance and/or the neutral gas velocity 

128 in the ionosphere-thermosphere system [Galand et al, 2001; Galand and Richmond, 

129 2001; Fang et al, 2007a, 2007b]. Both of these characteristics are crucial elements 

130 in the magnetosphere-ionosphere electrodynamics. Such wave-induced modification 

131 can be especially important equatorward of the low-latitude edges of the electron and 

132 proton auroral ovals where the wave-induced RC ion precipitation may be a dominant 

133 energy source. In addition, electrons and protons do not interact in the same way with 

134 the atmosphere. One should keep in mind that energetic protons ionize more efficiently 

135 than electrons do because their energy loss for each produced electron is smaller than 

136 that of energetic electrons [Galand et al, 1999]. Therefore, even if the proton energy 

137 flux is smaller compared to the electron flux, the response of the atmosphere to protons 

138 can be significant. The above arguments suggest that a self-consistent model of the 

139 magnetospheric electric field, RC, plasmasphere, and EMIC waves is needed to properly 



140 model wave spectral distribution and to improve the modeling of the large scale 
mi magnetosphere-ionosphere electrodynamics. 

142 In this study, we present a new computational model that is a result of coupling 

143 two RC models developed by our group. The first model deals with the large scale 

144 magnetosphere-ionosphere electro dynamic coupling and provides a self-consistent 
us description of RC ions and the magnetospheric electric field [ Liemohn et al, 2001 ; 
we Ridley and Liemohn, 2002 ; Liemohn et al, 2004 ], The second model is governed by a 
147 coupled system of the RC kinetic equation and the wave kinetic equation. This model 
ms self-consistently treats a mesoscale electrodynamic coupling of RC and EMIC waves, 

149 and determines the evolution of the EMIC wave power spectral density [ Khazanov et 

150 al, 2006 ; Khazanov et al., 2007 a]. The RC-EMIC wave model explicitly includes the 

151 wave growth/damping, propagation, refraction, reflection, and tunneling in a multi-ion 

152 magnetospheric plasma. Although RC ions and EMIC waves in the second model are 

153 treated self-consistently, the electric held is externally specified. So far, the above two 

154 models were used independently. As such, the main purpose of this paper is to present 

155 a new self-consistent model of the magnetospheric electric held, RC, plasmasphere, and 

1 56 EMIC waves along with initial results from the model simulations. The results presented 

157 in this study were obtained from simulations of the May 2 - 4 , 1998 geomagnetic storm, 

1 58 that we previously analyzed using an analytical formulation of the Volland-Stern electric 

159 held [ Khazanov et al., 2006 ; Khazanov et al, 2007 b]. 

leo This article is organized as follows: In section 2 we present a complete set 
i6i of the governing equations, and formulate the approaches used in the model 
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1 62 simulations. In the same section, we specify the initial/boundary conditions, and the 

163 interplanetary/geomagnetic characteristics, which drive our model. In section 3 the 

164 initial results from these simulations and discussion are provided. Finally, in section 4 

165 we summarize. 

lee 2. RC— EMIC Wave Model and Magnetosphere— Ionosphere 

167 Coupling 

lee 2.1. Governing Equations 

169 To simulate the RC dynamics we solve the bounce-averaged kinetic equation for 

170 the phase space distribution function of the major RC species (H + , 0 + , and He + ), 

171 as originally suggested in the models of Fok et al. [1993] and Jordanova et al. [1996]. 

172 The distribution function, F(r 0 , ip, E, p 0 ,t), depends on the radial distance in the 

173 magnetic equatorial plane r 0 , geomagnetic east longitude, kinetic energy E, cosine of 

174 the equatorial pitch angle /i 0 , and time t. For the F/e + -mode EMIC waves we also 

175 use the bounce-averaged kinetic equation. This equation describes a physical model of 

176 EMIC waves bouncing between the off-equatorial magnetic latitudes, which correspond 

177 to the bi-ion hybrid frequencies in conjugate hemispheres, along with tunneling across 

178 the reflection zones and subsequent strong absorption in the ionosphere (for the 

179 observational and theoretical justifications of this model see [ Gamayunov and Khazariov, 
iso 2008; Khazanov et al, 2007a]). The bounce-averaged wave kinetic equation was derived 
i8i in our previous paper [ Khazanov et al, 2006], and it explicitly includes the EMIC wave 
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182 growth/damping, propagation, refraction, reflection, and wave tunneling in a multi-ion 

183 magnetospheric plasma. In the present study, following Khazanov et al. [2006], we 

184 ignore the azimuthal and radial drifts of the wave packets during propagation, we do not 

185 include the wave tunneling across the stop zone, and consequently use a truncated wave 
lee kinetic equation. The resulting system of equations to drive RC-EMIC wave coupling 
i87 takes the form: 
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On the left-hand side of equation (1), all the bounce-averaged drift velocities are 
denoted as (• • •) and may be found in many previous studies [e. g., Khazanov et al, 
2003a]. The term on the right-hand side of this equation includes losses from charge 
exchange, Coulomb collisions, RC-EMIC wave scattering, and ion precipitation at low 
altitudes [e. g., Khazanov et al, 2003a]. Loss through the dayside magnetopause is 
taken into account, allowing a free outflow of the RC ions from the simulation domain. 
In equation (2), B w is the EMIC wave spectral magnetic field, u> and 9 0 are the wave 
frequency and equatorial wave normal angle, respectively, ( 9 0 ) is the bounce-averaged 
drift velocity of the wave normal angle, and (7) is a result of averaging the local 
growth/damping rate along the ray phase trajectory over the entire wave bounce period 
The factor (7) takes into account both the wave energy source due to interaction with 
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the RC ions and the energy sink due to absorption by thermal and hot plasmas. 

To perform bounce averaging in equation (2), the ray phase trajectory should 
be known, and we obtain it by solving the set of ray tracing equations. For a plane 
geometry these equations can be written as [e. g., Haselgrove, 1954; Haselgrove and 
Haselgrove, 1960; Kimura, 1966; Khazanov et al, 2006] 


dr ( dG/dk) r 

dt dG/du 

dX __{dG/dk) x 
dt dG/du 

dk r _ k d\ (dG/d r) r 
dt A dt dG/du 

dk\ k\dr {dG/d r) A 

dt r dt dG/du 


( 3 ) 

( 4 ) 

( 5 ) 

( 6 ) 


In equations (3)— (6) , the Earth-centered polar coordinate system is used to characterize 
any point P on the ray trajectory by length of the radius vector, r, and magnetic 
latitude, A. Two components, k r and k\, of the wave vector are given in a local Cartesian 
coordinate system centered on the current point P with its axes oriented along the 
radius vector and magnetic latitude direction, respectively. The function G (u, k, r) has 
roots for EMIC eigenmodes only, i. e., G — 0 at any point along the EMIC wave phase 
trajectories. Equations (3)-(6) are also used to obtain the off-equatorial power spectral 
density distribution for EMIC waves, which is needed to calculate the bounce-averaged 
pitch angle diffusion coefficient in the right-hand side of equation (1). (For more details 
about the system of equations (l)-(6) and its applicability please see our previous 
papers [Khazanov et al., 2003a; Khazanov et al, 2006; Khazanov et al, 2007a].) 


226 
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227 The bounce-averaged pitch angle diffusion coefficient on the right-hand side 

228 of equation (1) is a functional form of the EMIC wave power spectral density, and 

229 (7 (r 0 ,(p,t,u>,9 0 )) in equation (2) is a functional form of the phase space distribution 

230 function. So, there is a system of coupled equations, and the entire set of equations 

231 (l)-(6) self-consistently describes the interacting RC and EMIC waves in a quasilinear 

232 approximation. Compared to our previous RC-EMIC wave studies, which are based 

233 on equations (l)-(6) only [ Khazanov et al, 2006; 2007b], we are now going to take 

234 into account the magnetosphere-ionosphere coupling by self-consistently treating the 

235 current closure between RC and the ionosphere. 

236 Vasyliunas [1970] mathematically formulated a self-consistent model of the 

237 magnetosphere-ionosphere coupling by providing the basic equations governing the 

238 system. He outlined a logical chain of the model as follows: (1) the magnetospheric 

239 electric field determines the distribution of RC ions and electrons and, particularly, 

240 the total plasma pressure at any point; (2) from the plasma pressure gradients, the 

241 electric current perpendicular to the magnetic field can be calculated; (3) because the 

242 total current density should have zero divergence under magnetospheric conditions, the 

243 divergence of the perpendicular current density must be canceled by the divergence 

244 of FAC density, and so the divergence of the perpendicular current integrated along 

245 the entire field line gives the total FAC flowing into/ont of the conjugate ionospheres; 

246 (4) from the requirement that FAC is closed by the horizontal ohmic currents in the 

247 ionosphere, the distribution of the electric potential in the ionosphere can be found; 

248 and (5) the ionospheric potential can be mapped back into the magnetosphere along 
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geomagnetic field lines, and the requirement that this “new” magnetospheric electric 
held agrees with the “initial” magnetospheric held closes the magnetosphere-ionosphere 
system. 

To quantify the above logical chain, Vasyliunas [1970] used the following equations: 

j± ( r o> ¥>, S) = X (v J> ± + (B ■ V) B ) , (7) 

%(A(r„), 4 ?) = -B,(\(r„),ip) f VJ± — -da, (8) 

Jss B (r o, (p, s) 

VI,; = j ||, i sin x, I, = S x Bjj , (9) 

where P± and P\\ are the total plasma pressure (we neglect the electron pressure in the 
current study) perpendicular and parallel to the external magnetic held B, respectively, 
and J_|_ is the perpendicular current density. The FAC density at the ionospheric level 
is J\\ i (positive for current flowing into the ionosphere), Bi is the magnetic held in 
the ionosphere, and integration in equation (8) is done along the entire magnetic held 
line between foot points ss and sn . The coordinates (A (r 0 ) ,p) are the corresponding 
ionospheric latitude and MLT for the magnetic held line crossing the equatorial plane at 
(r 0 ,<p) (assuming that (p is the same at the equator and at the ionospheric altitude). In 
equations (9), I, and S are the height integrated horizontal ionospheric current density 
and conductivity tensor, respectively, and x is an inclination of the magnetic held (dip 
angle). The electric potential at the ionosphere level is <f>j, and V n is the velocity of the 
neutral gas in the ionosphere. Following many previous studies, in the present study we 
assume that the neutral gas corotates with the Earth and neglect the potential drop 
between the ionosphere and the equatorial magnetosphere [e. g., Ebihara et al, 2004], 
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272 Finally, it should be noted that, in general, equation (9) is written for the northern and 

273 southern ionospheres with the corresponding FAC juj, while equation (8) gives only 

274 the total FAC flowing into/out of the conjugate ionospheres but the obvious equation 


275 J||,j = j\\,i(ss) + j\\,i( s N ) is held. 

276 The set of equations (l)-(9) drives the RC, the EMIC waves, and the magnetospheric 

277 electric field in a self-consistent manner if all the initial and boundary conditions are 

278 specified and the ionospheric Hall and Pedersen conductances are known. A block 

279 diagram of the self-consistent coupling of the RC, EMIC waves, plasmasphere, and 

280 ionosphere is presented in Figure 1. The system characteristics in orange boxes are 

281 externally specified, and the dashed lines connect the model elements, which are 

282 currently not linked. 


Figure 1 


283 2.2. Approaches Used in Simulations 

284 The geomagnetic field used in the present study is taken to be a dipole field. It is 

285 a reasonable approximation for the present study because the most important results 

286 are obtained from simulations of the May 2-3, 1998 period (. Dst = -106 nT) when 

287 the Earth’s magnetic field is only slightly disturbed in the inner magnetosphere [e. g., 

288 Tsyganenko et al., 2003]. The convection electric field is calculated self-consistently as 

289 described in subsection 2.1, and the total electric field includes both the magnetospheric 

290 convection and corotation field. The equatorial cold electron density, n e , is obtained 

291 from the dynamic global core plasma model of Ober et al. [1997] . This model is basically 

292 the same as a time-dependent model of Rasmussen et al. [1993], which was used in our 
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293 previous studies, except the Ober et al. model is linked with a self-consistent electric 

294 held obtained from the system (1) — (9) , while the Rasmussen et ah model is driven by the 

295 Volland-Stern convection held [Volland, 1973; Stern, 1975] with Kp parameterization. 

296 Thus, the cold plasma density dynamics is also electrically self-consistent in our global 

297 RC-EMIC wave model. This is extremely important for a correct description of the 

298 EMIC wave generation/damping and propagation. In order to model the EMIC wave 

299 propagation and interaction with RC, we also need to know the density distribution 

300 in the meridional plane. In the present study we use a magnetic held model for the 

301 meridional density distribution, i. e., n e ~ B, because a more sophisticated analytical 

302 model by Angerami mid Thomas [1964] used in our previous studies [e. g., Khazanov 

303 et al., 2006] was found to give nearly the same results. The meridional model is 

304 then adjusted to the equatorial density model. So the resulting plasmaspheric model 

305 provides a 3D spatial distribution of the electron density. Besides electrons, the cold 

306 magnetospheric plasma is assumed to consist of 77% H + , 20% He + , and 3% 0 + , which 

307 are in the range of 10 — 30% for He + and 1 — 5% for 0 + following the observations by 

308 Young et al. [1977] and Horwitz et al. [1981]. Geocoronal neutral hydrogen number 

309 densities, needed to calculate loss due to charge exchange, are obtained from the 

310 spherically symmetric model of Chamberlain [1963] with its parameters given by Rairden 
3 n et al. [1986]. 

312 During the main phase of major storms, RC 0 + may dominate [e. g., Hamilton et 

313 al., 1988; Daglis, 1997] and, as a result, contribute to strong damping of the 77e + -mode 

314 EMIC waves [ Thorne and Horne, 1997]. Although there is no doubt that, in principle, 
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315 this process is important, let us evaluate the validity of excluding the He + - mode 
3 ie damping by RC 0 + in the May 2-4, 1998 storm simulation. Using the RC kinetic model 

317 of Jordanova et al. [1998], Farrugia et al. [2003] found that during the main phase 

3 1 8 of the May 4, 1998 storm the energy density of RC H + is greater than twice that of 

3 1 9 0 + at all MLTs, and the contribution of He + to the RC energy content is negligible. 

320 This implies that the RC 0 + content does not exceed 30% during the main phase of 

321 this storm. This estimate was obtained from a global simulation, which did not include 

322 oxygen band waves. On the other hand, Brdysy et al. [1998] observed a very asymmetric 

323 0 + RC during the main phase of the April 2-8, 1993 storm, which may suggest that a 

324 majority of the RC oxygen ions get lost before they reach the dusk MLT sector. This 

325 result is difficult to explain in terms of charge exchange and Coulomb scattering, and 

326 suggests that the production of EMIC waves contributes significantly to RC 0 + decay 

327 during the main and early recovery phases. In other words, due to the generation of 

328 the 0 + -mode EMIC waves, most RC 0 + might precipitate before reaching the dusk 

329 MLT sector [Brdysy et al., 1998]. Therefore, to estimate the RC 0 + content correctly, 

330 the 0 + -mode should be included in the simulation, and it is likely that Farrugia et al. 

331 [2003] overestimated the RC 0 + content during May 4, 1998. Moreover, the calculations 

332 of Thorne and Horne [1997] clearly demonstrated that even the RC 0 + percentage 

333 noted above cannot significantly suppress the L/e + -mode amplification, and only slightly 

334 influences the resulting growth; inclusion of 26% 0 + in the RC population causes the 

335 net wave gain to decrease by only 20%. In addition, the most important results shown 

336 in the present study are obtained from simulations of the May 2-3, 1998 period, i. e., 
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337 the first main (. Dst = -106 nT) and recovery phases of the May 1998 large storm, when 

338 the RC 0 + content should be even smaller than the Farrugia et ah estimate for May 4 , 

339 1998 . It is for these reasons that we chose to exclude RC 0 + in the present simulations, 

340 and to assume that the RC is entirely comprised of energetic protons. 

341 Equation ( 9 ) must be solved taking into account the contributions from both 

342 the northern and southern ionosphere. Because in the present study we assume the 

343 magnetic field lines to be equipotentials, the northern and southern ionospheres can 

344 just be replaced by an effective single ionosphere with E = E^ + E N , and total FAC 

345 flowing into/out of it. After the resulting equation is solved, and is found, we 

346 can easily calculate the FACs j\\ t i(ss) and j\\ t i(sN) flowing into/out the southern and 

347 northern ionosphere. 

348 The ionospheric Hall and Pedersen conductances in our model are not calculated 

349 self-consistently but rather specified by empirical models. The resulting conductance 

350 arises from four sources: ( 1 ) direct solar extreme ultraviolet (EUV), ( 2 ) scattered solar 

351 EUV on both sides of the terminator, ( 3 ) starlight, and ( 4 ) auroral particle precipitation. 

352 The direct solar conductance is controlled by the solar zenith angle and the solar UV 

353 and EUV radiations, which correlate with the solar radio flux index F10.7. In the present 

354 study we use the empirical model of Moen and Brekke [ 1993 ] for determining direct 

355 solar conductance. The scattered solar EUV and starlight conductance models are taken 

356 from the study of Rasmussen and Schunk [ 1987 ]. In order to specify the conductance 

357 from auroral precipitation, we use either the Hardy et al. [ 1987 ] statistical model or an 

358 empirical relationship between the FACs and the local Hall and Pedersen conductance 
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359 established by Ridley et al. [2001; 2004], The Hardy et al. model is compiled from 

360 the electron precipitation patterns obtained by the DMSP satellites and gives the Hall 

361 and Pedersen conductance as a function of MLT and magnetic latitude for seven levels 

362 of activity as measured by Kp. The Ridley et al. relationship was derived using the 

363 assimilative mapping of ionospheric electrodynamics (AMIE) technique [ Richmond 

364 and Kamide, 1988]. The AMIE technique was run at a one-minute cadence for the 

365 entire month of January 1997, using 154 magnetometers. This resulted in almost 

366 45000 2D maps of the Hall and Pedersen conductances and FAC. The conductance was 

367 derived from the Ahn et al. [1998] formulation, which relates ground-based magnetic 

368 perturbations to the Hall and Pedersen conductances. The Ridley et al. analysis showed 

369 an exponential relationship between the local FAC and the conductance [see Amm, 

370 1996; Goodman, 1995]: 

371 E = Eoe^l'M, (10) 

372 where the constants E 0 and A are independent of the magnitude of jy*, but depend 

373 on location and whether the current is upward or downward. Although the Ridley 

374 et al. relationship is entirely empirical and not based on Erst principles, by using it 

375 we introduce into the model at a degree of self-consistency between the ionospheric 

376 conductance and FAC. This is a principle modification because a self-consistent 

377 description of the ionospheric conductance makes equation (9) nonlinear compared 

378 to the case of statistical conductance model. For previous use of the Ridley et al. 

379 relationship in the RC simulation see Liemohn et al. [2005]. 
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380 To conclude this subsection, we note that the numerical implementations used to 

381 solve equations (l)-(6) are described in details in our previous publications [ Khazanov 

382 et al, 2003a; 2006], and to solve equation (9) a preconditioned gradient reduction 

383 resolution (GMRES) solver is used [Ridley et al, 2004], The GMRES method is robust 

384 enough to handle a wide variety of FAC and conductance patterns. 

385 2.3. Initial and Boundary Conditions 

386 The initial RC distribution is constructed from the statistically derived quiet time 

387 RC proton energy distribution of Sheldon and Hamilton [1993] and the pitch angle 

388 characteristics of Garcia and Spjeldvik [1985]. The night-side boundary condition 

389 for equation (1) is imposed at the geostationary distance, and it is obtained using 

390 flux measurements from the Magnetospheric Plasma Analyzer [ Bame et al, 1993] 

391 and the Synchronous Orbit Particle Analyzer [ Belian et al, 1992] instruments on the 

392 geosynchronous LANL satellites during the modeled event. Then, according to Young et 

393 al [1982] and Liemohn et al. [1999], we divide the total flux measured at geostationary 

394 orbit between the RC H + , 0 + , and He + depending on geomagnetic and solar activity 

395 as measured by Kp and T10.7 indices. Only the H + flux is used as a boundary condition 

396 in the simulation. 

397 In the present study, the poleward boundary for equation (9) is taken at magnetic 

398 latitude A = 69°. On this boundary, we specify the electric potential using either 

399 the Weimer [1996] statistical model (hereinafter the W96 model), which is driven by 

400 the interplanetary magnetic field (IMF) By, Bz components and solar wind velocity, 
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or the convection model of Volland and Stern [Volland, 1973; Stern, 1975] with 
Kp parameterization given by Maynard and Chen [1975] and shielding factor of 2 
(hereinafter the VS model). The second boundary condition is specified at A = 30°, and 
we use either the W96 model or the VS model, both of which give the potential close 
to zero at that latitude. It should be noted that the result of calculation is insensitive 
to the choice of the lower boundary condition, as demonstrated by Wolf [1970]. So, the 
magnetospheric electric held is calculated self-consistent ly in the domain 30° < A < 69°. 
At the same time, we should emphasize that, compared to RC, the cold electron density 
is modeled in a more extended domain of L < 10, and in order to specify the electric 
held in the entire L < 10 region, we use either the W96 or the VS model for the 
magnetic latitude above A = 69°. 

The initial RC, plasmasphere, and EMIC wave distributions are derived 
independently and, moreover, they have nothing to do with a particular state of the 
magnetosphere/plasmasphere system during a simulated event. Only the boundary 
conditions provided by the LANL satellites can be considered as data reflecting a 
particular geomagnetic situation (and, to a certain extent, the employed ionospheric 
conductance model and an imposed cross polar cap potential drop). Therefore, before 
the simulation of a particular geomagnetic event can occur, we first must find an 
appropriate initial state for the RC, electric held, plasmasphere, and EMIC waves 
that is self-consistent and reflects the particular geomagnetic situation. To obtain 
the self-consistent initial distributions for the entire system, we first prepared the 
plasmasphere by running the Ober model for 20 quiet days. Then, at 0000 UT on 
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1 May, 1998, a simulation of equations (1)-(10) was started using all the controlling 
parameters and the initial/boundary conditions along with a background noise level for 
the He + - mode EMIC waves [e. g., Akhiezer et al, 1975]. We ran the model code for 
24 hours to achieve a quasi-self-consistent state for the system. Note that 24 hours 
has nothing to do with the typical time for wave amplification and instead reflects 
the minimum time needed to adjust the RC and waves to each other and to the real 
prehistory of a storm. The self-consistent modeling of the May 1998 storm period 
was started at 0000 UT on 2 May (24 hours after 1 May 0000 UT) using solutions of 
equations (1), (2), and the cold plasma distribution at 2400 UT on 1 May as the initial 
conditions for further simulation. 


2.4. Interplanetary and Geomagnetic Drivers for the Model 


The ionospheric boundary condition in our simulations is driven either by IMF By , 
Bz components and solar wind velocity (the W96 model) or the 3-hour Kp index (the 
VS model). The Hardy et al. [1987] ionospheric conductance model is driven by Kp. 
All of these driving parameters are shown in Figure 2 during the May 2-4, 1998 period. 
Interplanetary data are obtained from the Magnetic Field Investigation [ Lepping et al, 
1995] and the Solar Wind Experiment [ Ogilvie et al, 1995] instruments aboard the 
WIND satellite. The interplanetary configuration of May 1-5, 1998 consists of a coronal 
mass ejection (CME) interacting with a trailing faster stream [ Farrugia et al, 2003]. 
The CME drives an interplanetary shock observed by the instruments aboard the WIND 
spacecraft at about 2220 UT on May 1. Three episodes of the large negative IMF Bz 
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444 component were monitored. The first episode started at ~ 0330 UT on May 2 (27.5 

445 hours after May 1, 0000 UT), the second at 0230 UT on May 4 (74.5 hours after May 

446 1, 0000 UT), and the third (not shown) at ~ 0200 UT on May 5 (98 hours after May 

447 1, 0000 UT). These caused a “triple-dip” storm with the minimums Dst = —106 nT, 

448 Dst = —272 nT, and Dst = —153 nT (not shown). The planetary Kp index reached 

449 maximum values of Kp « 7~ and Kp 9~ at the times when Dst minimums were 

450 recorded. 

451 3. Results and Discussion 

452 3 . 1 . Magnetospheric Electric Field 

453 The cross polar cap potential (CPCP) drop gives a rough quantitative assessment 

454 of the strength of convection in the inner magnetosphere. We calculate the CPCP drop 

455 as a difference between the maximum and minimum values of the potential at A = 67.5° 

456 (at L « 7). Results of our calculations are shown in Figure 3. The lines in red, green, 

457 and blue show results from a self-consistent simulation, while the CPCP drop shown in 

458 black is for reference purposes only. Note that the red line lies somewhat higher than 

459 the black one. This is because we do not calculate FACs between A = 69° and A = 67.5° 

460 in the present simulations, and so there is no shielding taken into account unlike in 

461 the analytical formulation of the VS potential (black line in Figure 3). When the W96 

462 model is imposed at A = 69°, the CPCP drops are very similar for both conductivity 

463 models, and the blue line is just slightly higher than the green one. The CPCP drop 
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464 resulting from the VS model is larger during the majority of May 2-4, except for about 

465 13 hours on May 2 and 12 hours on May 4, when the CPCP drop from the W96 model 

466 is greater. It is seen that the W96 potential drop spikes to 300 kV during the main 

467 phase on May 4, whereas the VS boundary condition results in a maximum CPCP drop 

468 of only 150 kV. 

469 Although the CPCP drop may serve as an overall measure of the convective 

470 strength, it does not give the morphology and strength of the electric held in the inner 

471 magnetosphere. To provide such insight, we selected six snapshots of the equatorial 

472 electric held patterns from May 2, and one snapshot at hour 77 (0500 UT on May 4). 

473 The corresponding electric potential contours are shown in Figure 4. The view is over 

474 the North Pole with local noon to the left. We present results for three runs. The 

475 equipotentials from a simulation with the VS model at the high latitude ionospheric 

476 boundary and the Hardy et al. conductance are shown in the hrst row. The other 

477 two runs are performed with the W96 model applied at A = 69°, and differ only by 

478 the conductance model assumed. The second row shows results for the Hardy et al. 

479 conductance model, while the third row is for a case when the Ridley et al. empirical 

480 relationship between the FAC and conductance is used. The potential configurations 

481 in Figure 4 are similar to those from the Rice Convection Model [e. g., Garner et al, 

482 2004], Overall, there are qualitatively the same large-scale potential distributions in 

483 all three models, presented in Figure 4 with a well defined large-scale dawn-to-dusk 

484 electric held. Despite this, the potential patterns reveal large differences in both the 

485 magnitude of the potential and the shape of the contours. This suggests a difference in 
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486 the fine structure of the electric field distribution since this field is proportional to the 

487 gradient of the potential. 

488 One obvious feature observed in Figure 4 is a significantly enhanced electric field 

489 in the region L « 3 — 4 in the dusk-post-midnight MLT sector at hour 77 (and, not 

490 shown, at hour 76). This radially narrow intensification of the radial electric field 

49 1 (poleward electric field in the ionosphere) creates a westward flow channel, mainly in 

492 the dusk-to-midnight MLT sector, while a region of westward (antisunward) convection 

493 is also observed in the post-midnight sector equatorward of L = 3 (see Figure 4). This 

494 westward flow channel has come to be called the subauroral polarization stream (SAPS) 

495 [ Foster and Burke, 2002; Foster and Vo, 2002], The SAPS effect arises from the region 2 

496 FACs, which flow down into the subauroral ionosphere and close the region 1 FACs 

497 through the poleward Pedersen currents. Because of the low conductance at subauroral 

498 latitudes, the Pedersen current generates an intense poleward electric field between the 

499 region 2 FAC and the low-latitude edge of the auroral particle precipitation [Southwood 

500 and Wolf, 1978; Anderson et al, 1991, 1993; Ridley and Liemohn, 2002; Mishin and 

501 Burke, 2005], 

502 To show the potential structure and electric field inside the SAPS region, we took 

503 two meridional cuts across the entire simulation domain and the corresponding results 

504 are shown in Figure 5. Figures 5a, b show the potential profiles on the dawn-dusk 

505 meridian for hours 33 and 77. Results for three simulations are presented along with a 

506 profile for the analytical VS model. The corresponding equatorial radial electric fields 

507 are shown in Figures 5c, d for MLT=18. Only a slight electric field intensification 
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508 (< 2.7 mV/m) is observed in the dusk sector for hour 33 (see Figure 5c), while we see 

509 an extremely developed SAPS in Figure 5d (< 13.4 mV/m). The strongest electric 

510 held intensification in Figure 5d takes place for cases when the W96 model is used 

5 11 in combination with either the Hardy et al. conductance model or the Ridley et al. 

512 relationship. In the latter case, we see a slightly stronger electric held in the dusk MLT 

513 sector and a developed dawnside electric held of about 5 mV/m (see Figure 5b). 

514 Although the SAPS localization is correctly predicted by our model, it is likely that 

515 the SAPS electric held in Figure 5d is overestimated for the W96 boundary condition. 
5 ie Indeed, from the statistical model based on the electric held data measured by the 

517 Akebono/EFD instrument, Nishimura et al. [2007] derived the equatorial Ey electric 

5 1 8 held component in the dusk SAPS region to be 6 mV/m during the main phase of storm. 

519 It should be noted, however, that the SAPS electric held can sometimes reach more 

520 than 10 mV/m during the main phase of geomagnetic storms [Shinbori et al, 2004], and 

521 the CPCP drop derived by Nishimura et al. [2007] is 180 kV, whereas in our simulation 

522 it is 300 kV. The measurements taken by the double-probe electric held instrument 

523 on-board the CRRES spacecraft show a similar electric held magnitude [Wygant et al, 

524 1998]. There are at least two reasons that may lead to an overestimation of the SAPS 

525 electric held in our simulations. (1) Because the W96 model was constructed from data 

526 with IMF under 10 nT, this model essentially overestimates the CPCP drop during the 

527 May 4 event when IMF was around 40 nT [e. g., Burke et al, 1998]. (2) In the present 

528 simulations, we did not take into account the FACs beyond geostationary orbit, which 

529 may contribute essentially to the shielding of midlatitudes from a high latitude driving 
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530 convection field; the effect of FAC is proportional to the volume of the magnetic flux 

531 tube, and from the estimate by Vasyliunas [1972] the effect of FAC at L=6.6 is about 

532 20% of the FAC effect at L=10. Both of these issues will be addressed in future studies. 

533 3 . 2 . Plasmasphere 

534 The plasmapause, and/or dayside plume, and/or detached plasma are the favorable 

535 regions for EMIC wave generation in the inner magnetosphere. This is because 

536 the density gradient there is enhanced and counteracts refraction caused by the 

537 magnetic field gradient and curvature [e. g., Horne and Thorne, 1993; Fraser et al, 

538 2005; Khazanov et al., 2006]. As a result, the net refraction is suppressed at the 

539 plasmapause/plumc edge allowing wave packets to spend more time in the phase region 

540 of amplification. Thus, the cold plasma distribution is extremely crucial for EMIC 

541 wave excitation. Both the convection and the corotation electric fields control the cold 

542 plasma dynamics. As such, we will first present the snapshots of the total electric 

543 potential obtained from our simulations. Figure 6 shows the resulting equipotential 

544 contours, that also coincide with the instantaneous cold plasma flow. The most striking 

545 reconfiguration of the potential is observed in the second and third rows in the 28 and 

546 30 hour snapshots. Referring to Figure 3, we see that starting at hour 28 the CPCP 

547 drop increases by about 100 kV during one hour for the W96 convection model. The 

548 strong convection causes a shrinking of the closed equipotential contours as shown in 

549 Figure 6 (there is stronger shrinking during hour 29). Later, an extremely developed 

550 SAPS is observed at hours 76-77 (see subsection 3.1), and the overshielding electric field 
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551 (negative E Y ) following a decrease of the CPCP difference in the W96 model is found 

552 in the inner magnetosphere at hour 79 (not shown). 

553 Figure 7 shows the selected distributions of the equatorial cold plasma density for 

554 three self-consistent simulations. For each run, the plasmasphere was first prepared by 

555 running the Ober code for 20 quiet days. Then, starting at 0000 UT on 1 May, 1998, 

556 we solved the equations (1)-(10) using the initial and boundary conditions and the 

557 time series for all controlling parameters (see subsection 2.3). For the VS model (first 

558 row), a broad dayside plume is formed a few hours before hour 28. Subsequently, up to 

559 hour 39 gradual intensification of the convection (see Figures 3 and 4) causes nightside 

560 plasmaspheric erosion and the plume narrowing in the MLT extent. The latter takes 

561 place mostly in the eastward flank of the plume where the convection and corotation 

562 fields reinforce each other, while the duskside plume edge remains roughly stationary 

563 [Spasojevic et al, 2003; Goldstein et al, 2005]. During the following storm progression, 

564 the magnetospheric convection held driven by the VS potential drop remains relatively 

565 high (see Figure 3), and the convection patterns are relatively steady (3-hour cadence). 

566 Compared to the second and third rows in Figure 7, these result in the most eroded and 

567 shrunken plasmasphere at hour 77 with a well-defined nightside plasmapause (compare 
see these results with Figure 7 in [Khazanov et al, 2006] where the entire plasmasphere was 

569 driven by the analytical formulation of the VS potential). 

570 Cold plasma density distributions in the second and third rows of Figure 7 are 

571 qualitatively similar to each other, but exhibit quite a bit of difference compared to 

572 distributions in the first row. At hour 28, the plasmasphere is well-populated, and the 
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573 plasmapause is well-defined. Starting at hour 28, an increase of the CPCP drop by 

574 100 kV during one hour (see Figure 3) causes formation of the plume by hour 29 (not 

575 shown), and the presented snapshots at hour 30 are close to those at hour 29. One of the 

576 most distinguishable features observed in the second and third rows is the presence of a 

577 cold plasma on the nightside. To emphasize the existence of the recirculated detached 

578 plasma material, we show in Figure 8 the detailed plasma density evolution in the 

579 extended domain of L < 10. It is clearly seen in Figure 8 how this recirculated detached 

580 plasma is forming and reentering the inner magnetosphere. The radial electric field for 

581 MLT=18 and 19 is also shown in Figure 9 for hours 28 and 29. The negative electric field 

582 in the outer region in Figure 9b is resulting in plasma recirculation. However, we have 

583 to emphasize that a great care is needed to interpret these simulation results. During 

584 an extreme condition, the W96 model may predict a two-cell convection pattern with 

585 its focuses located at low latitude. The anti-sunward ionospheric plasma flow predicted 

586 by the W96 model may correspond to the lobe and the outer part of low-latitude 

587 boundary layer (LLBL) in the magnetosphere. In the daysicle magnetosphere, when the 

588 plasmaspheric cold plasma is transported to LLBL, the cold plasma will flow in the 

589 anti-sunward direction [e. g., Ober et al, 1998]. At the same time, reentry of the cold 

590 plasma from LLBL back to the magnetosphere may not be simple as predicted by the 

591 W96 model. 

592 Although the cold plasma recirculation is seen in both the second and the third 

593 rows of Figure 7, the observed similarity is only qualitative and all the quantitative 

594 characteristics are quite different. After hour 39, the W96 CPCP drop decreased and 
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595 fluctuated around 50 kV except for four hours on May 4 when the CPCP drop spikes 

596 to 300 kV during the second main phase of the storm (see Figure 3). In both cases, 

597 the resulting plasmaspheres at hour 77 are extremely diffusive with shallow density 

598 gradients. This is because the anti-sunward plasma flow is especially strong during the 

599 second main phase of the storm. To demonstrate that, we show in Figure 10 the total 

600 radial electric field versus MLT for L=8, 9, and 10 at hour 77. The negative radial 

601 electric field in the afternoon-premidnight MLT sector causes a counter clockwise plasma 

602 convection. The MLT extent of the negative electric field in the afternoon-premidnight 

603 MLT sector grows with L-shcll, resulting in the backward plasma flow for MLT > 15 

604 at L=10. This recirculation supplies the cold plasma in the nightside preventing the 

605 plasmasphere to be eroded. At the same time, as we emphasized above, a great care is 
6oe needed to interpret these results. 

607 To show the equatorial cold plasma density profiles during the periods of a 
60s well-defined and a shallow plasmapause we selected hours 33 and 77. Results of our 

609 simulations are shown in Figure 11. We see a “classical” profile of the plasmapause 

610 for hour 33, when the plasma density decreases about two orders of magnitude over 
6n 0.5 — 0.75 Re- The combination of the W96 model and the Ridley et al. relationship 

612 results in a detached plasma with a peak density of 20 cm -3 , which is clearly observed 

613 in Figure 11a (see also the third row in Figure 7). During hour 77, the plasmasphere 

614 driven by the VS CPCP drop is the most eroded and, although the plasmasphere 

615 boundary layer is wider than in Figure 11a and the plasma density drop is smaller, 

616 the plasmapause is still well-defined. For simulations with the W96 potential at the 
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617 high latitude ionospheric boundary, both density profiles shown in Figure lib exhibit a 

618 shallow density gradient without the plasmapause while there is a clear change of the 

619 profile slope for the W96-Hardy et al. result. Note that there are also no steep density 

620 gradients outside of geostationary orbit (not shown). 


62 i 3.3. RC Proton Precipitation 


622 The convection electric field controls the global precipitating patterns of RC. As 

623 RC protons approach the Earth via the convection electric field, they precipitate into 

624 the loss cone because the equatorial loss cone angle increases with decreasing L-shell 

625 somewhat more than the equatorial pitch angle increases [e. g., Jordanova et al, 1996]. 

626 Note that precipitation due to Coulomb collisions with thermal plasma takes place 

627 mainly inside the plasmapause, and the wave-induced ion precipitation is organized in 

628 the radially narrow regions in the plasmasphere boundary layer [e. g., Gurgiolo et al., 

629 2005; Khazanov et al, 2007b]. The RC proton precipitating fluxes integrated over two 

630 energy ranges 1 — 50 keV and 50 — 400 keV are calculated as 

1 rE 2 r 1 r 1 

631 J ic = — / &E / d/ioj, flic = / d/i 0 , (11) 

' he J El J flic ^ flic 


632 where ni c is the cosine of the equatorial pitch angle at the boundary of the loss cone, 

633 and j is the equatorial differential flux of RC protons. The snapshots of the fluxes for 

634 low and high energies are shown in Figure 12 and 13, respectively. The results from 

635 three self-consistent runs with a specified combination of the high latitude ionospheric 
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boundary potential and conductance model are shown. For low energy, the most intense 
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637 precipitating fluxes near the end of the second main phase (hour 77) are observed in 

638 the second and third rows of Figure 12 when the W96 model is used. This takes place 

639 because the convection held is strongest in these two cases (see Figure 4). The spot-like 

640 spatial structure in the postnoon-midnight MLT sector is due to the wave-induced 

641 precipitation with the strongest fluxes up to 10' cm _2 s _1 sr _1 . 

642 The penetrating electric held driven by the W96 boundary held causes precipitation 

643 of energetic RC ions well earthward of the low energy ion precipitation. It is clearly 

644 seen in Figure 13 that the W96 boundary potential leads to a strong precipitation of 

645 the high energy ions near the inner edge of RC during the second main phase on May 4. 

646 The high energy precipitating fluxes maximize at about two times stronger magnitude 

647 than the maximal fluxes observed in the range 1 — 50 keV. 


648 3.4. Energy Distribution for He + — Mode EMIC Waves 

The coupling of the magnetosphere and ionosphere by the region 2 FACs gives a 
self-consistent description of the magnetospheric electric held. This held controls the 
convective patterns of both RC ions and the cold plasmaspheric plasma, changing the 
conditions for EMIC wave generation/amplihcation. The equatorial (MLT, L-shell) 
distribution of the squared wave magnetic held, 


CWmax r 7T 

By, (r o, = / du d6 0 B^ (r 0 , </?, t, u, 8 0 ) , 

J ^min 

649 is shown in Figure 14 for the He + - mode EMIC waves. As before, the results from three 

650 self-consistent simulations are presented. Comparing Figure 14 with the cold plasma 
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651 density distribution in Figure 7, we see that EMIC waves are distributed in the narrow 

652 regions inside the plasmasphere boundary layer where the density gradient is enhanced. 

653 Although, during hours 30-39, the spatial wave distributions in the first and second 

654 rows look similar, on average, there are much more waves in a simulation with the VS 

655 boundary condition than in a simulation with the W96 potential during entire May 2. 

656 Moreover, there are practically no waves in the latter simulation after hour 39 (not 

657 shown) while in the former case we observe the extended regions of intense waves during 

658 the majority of the time up to hour 60 (not shown). This is because the plasmapause 

659 is well-defined and the CPCP drop is higher in the case of the VS potential boundary 

660 compared to the case of the W96 potential when the plasmasphere is highly diffusive (a 

661 shallow density gradient) and RC is less intense (lower the local growth rate). 

662 The density distributions in the second and third rows of Figure 7 demonstrate quite 

663 a bit of difference in the after-dusk MLT sector starting at hour 33. The plasmapause 

664 in the third row is located closer to the Earth, and the density gradient is shallowed 
ees by the detached plasma. At the same time, we observe much less wave activity in the 
eee third row of Figure 14 than in the second row. This is likely due to the effect of the 
667 density distribution, because the global potential drop is even higher in the third row of 
ees Figure 4 (suggesting a more intense RC) compared to the second row. 

669 There are practically no waves during the second main and recovery phases, 

670 except for moderate wave activity in the hour 77 snapshots in the Erst and third rows 

671 of Figure 14. In the case of the VS-Hardy et al. combination, the plasmapause is 

672 well-defined during hour 77 (see Figures 7 and 11) and waves can grow despite a less 
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673 intense RC in this case. On the other hand, the RC is strongly developed for the case 

674 of the W96 potential, and wave growth rate is essentially higher than in the first row, 

675 causing a wave generation despite the plasmasphere being extremely diffusive and the 

676 density gradient being shallow. 

677 3.5. Ionosphere Reconfiguration and Communication Time 

678 All of the results presented above were obtained from simulations when only a 

679 30 min time delay between WIND and the high latitude ionospheric boundary was 

680 applied. Both the reconfiguration time needed to reestablish a new potential pattern 

681 throughout the ionosphere and communication time between the ionosphere and the 

682 equatorial magnetosphere were assumed to be zero. These allowed us to update the 

683 equatorial electric field for each time step (a minute). However, this is not the case 

684 and both the ionospheric reconfiguration time and the Alfven propagation time are 

685 essentially higher than a minute [e. g., Ridley et al., 1998]. This implies that the 
ese ionosphere cannot reconfigure instantly in response to change of the interplanetary 
687 conditions, and that the magnetospheric electric field requires a finite time to be 
ess reestablished. 

689 Ridley et al. [1998] studied the ionospheric convection changes associated with 

690 changes of the IMF. They found that the total reconfiguration time of the ionosphere is in 

691 the range 3-26 min with an average of 13 min. Taking 7 min as a typical communication 

692 time between the ionosphere and the equatorial magnetosphere (for example, the 

693 magnetopause-ionosphere communication time is 8.4 ± 8.2 min as estimated by Ridley 
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694 et al. [ 1998 ]), on average, the same 13 min are needed to reestablish a new potential 
ess pattern in the magnetosphere but a 7 min delay should be applied to the ionospheric 

696 pattern. Because a great deal of scatter was reported for both time scales, below we 

697 simply adopt 20 (= 13 + 7 ) min as a time needed to reestablish a new potential pattern 

698 in the equatorial magnetosphere. 

699 To assess the importance of the finite ionospheric reconfiguration and communication 

700 time effect, we reran the “W 96 -Hardy et al.” simulation. Starting at hour 24 , we 

701 averaged the interplanetary parameters and FACs over a 20 min window before passing 

702 them to the ionospheric solver, and updated the equatorial electric field only once every 

703 20 min. Figure 15 shows the equatorial potential contours from this simulation along 

704 with the contours from the previous simulation, when the equatorial electric field is 

705 updated for each time step. The results during seven consecutive hours are shown 
7oe (hours 35 - 41 ). The potential distributions in the first and second rows are quite a 
707 bit different suggesting that the finite ionospheric reconfiguration and communication 
70 s time effect may be important, especially for the fine temporal-spatial structure of 

709 the plasmasphere-magnetosphere system. Although the “new” electric field alters the 

710 RC, wave, and cold plasma distributions, we show only the results for cold plasma 
7n density. Figure 16 demonstrates a difference in the cold plasma density distribution 

712 introduced by the effect of a finite time required to reestablish a “new” distribution 

713 of the magnetospheric electric held. Although the density distributions in these two 

714 simulations are identical at hour 24 , the plasmapause/plume shapes get a visible 

715 difference in the dawn-noon MLT sector starting at hour 29 (not shown). Later, starting 
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7ie at hour 35 , an essential difference between the density distributions is observed in the 

717 night MLT sector (see Figure 16 ). After hour 56 , the cold plasma density distributions 

718 in these two simulations are similar. This is expected after a longterm interval of system 

719 evolution, while the fine density structure still differs from time to time depending on 

720 the differences in the electric field distributions in these two simulations. 

721 Although a more sophisticated methodology is required to treat and separate the 

722 effects of the finite ionospheric reconfiguration and communication time, Figures 15 and 

723 16 clearly demonstrate that the finite time effect is important, especially for the fine 

724 temporal-spatial structure of the system. This implies that the instant interplanetary 

725 parameters cannot be used in order to specify the outer ionospheric boundary condition, 

726 but rather some kind of the averaging procedure should be applied to these parameters 

727 before passing them to the ionospheric solver. 

728 4. Summary 

729 The scattering rate of magnetospheric RC ions and relativistic electrons by EMIC 

730 waves is not only controlled by the wave intensity-spatial-temporal distribution but 

731 strongly depends on the spectral distribution of the wave power. There is growing 

732 experimental [Anderson et al, 1996 ; Denton et al, 1996 ; Anderson et al, 1992 b; Fraser 

733 and Nguyen, 2001 ; Meredith et al., 2003 ] and theoretical [ Horne and Thorne, 1993 ; 

734 Khazanov et al, 2006 ] evidence that EMIC waves can be highly oblique in the Earth’s 

735 magnetosphere. Compared to field-aligned waves, the highly oblique wave normal 

736 angle distributions can dramatically change the effectiveness (an order of magnitude 
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or more) of both the RC proton precipitation [Khazanov et al., 2007b] and relativistic 
electron scattering [ Glauert and Horne, 2005; Khazanov and Gamayunov, 2007]. 
Strong sensitivity of the scattering rates to the wave spectral characteristics suggests 
that in any effort to model EMIC wave-induced heating and/or scattering of the 
magnetospheric particles, the wave spectral distribution requires special care and should 
be properly established. Unfortunately, there are still very few satellite-based studies 
of EMIC waves, especially during the main phase of magnetic storms, and currently 
available observational information regarding EMIC wave power spectral density is poor 
[Engebretson et al, 2008]. So, a combination of comprehensive theoretical models and 
available data should be utilized to obtain the power spectral density of EMIC waves 
on the global magnetospheric scale throughout the different storm phases. To the best 
of our knowledge, there is only one model that is able to simulate a spatial, temporal 
and spectral distribution of EMIC waves on the global magnetospheric scale during the 
different storm phases [ Khazanov et al., 2006]. This model is based on first principles 
and is governed by a coupled system of the RC kinetic equation and the wave kinetic 
equation, explicitly including the wave generation/damping, propagation, refraction, 
reflection and tunneling in a multi-ion magnetospheric plasma. 

The convective patterns of both the RC ions and the cold plasmaspheric plasma 
are controlled by the magnetospheric electric field, thereby determining the conditions 
for interaction of RC ions and EMIC waves. Therefore, this electric field is one of 
the most crucial elements in simulating the wave power spectral density on a global 
magnetospheric scale. Self-consistent simulation of the magnetosphere-ionosphere 
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759 system should provide, at least in principle, the most accurate theoretical electric 

760 field [Vasyliunas, 1970; Jaggi and Wolf, 1973]. The need for a self-consistent model 

761 of the magnetospheric electric field, RC, plasmasphere, and EMIC waves is evident. 

762 In the present study we have incorporated the large scale magnetosphere-ionosphere 

763 electrodynamic coupling in our previous self-consistent model of interacting RC ions 

764 and EMIC waves [Khazanov et al., 2006]. The resulting computational model treats 

765 self-consistently not only EMIC waves and RC ions but also the magnetospheric electric 

766 field, RC, and plasmasphere. 

767 A few runs of this new model were performed to get a qualitative assessment of 

768 the effects of the high latitude ionospheric boundary condition and the ionospheric 

769 conductance. The results presented in this study were obtained from simulations 

770 of the May 2-4, 1998 geomagnetic storm (mostly the May 2-3 period). We have 

771 performed three simulations that differ by the electric potential specified at the high 

772 latitude ionospheric boundary (we used the W96 model and the VS model with Kp 

773 parameterization), and/or the ionospheric conductance from auroral precipitation 

774 (utilizing the Hardy et al. conductance model and the Ridley et al. relationship between 

775 the FACs and the conductance). The following three combinations have been used in 

776 the simulations: (1) the VS model and the Hardy et ah model; (2) the W96 model and 

777 the Hardy et ah model; and (3) the W96 model and the Ridley et ah relationship. In 

778 addition, one more simulation has been done: (4) the W96 model and the Hardy et 

779 ah model applying a 20 min window as the time needed to reestablish a new potential 

780 pattern in the magnetosphere. The RC in the present study has been simulated inside 
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781 geostationary orbit only, and the high latitude ionospheric boundary has been placed 

782 near the ionospheric projection of this orbit. The findings from our initial consideration 

783 can be summarized as follows: 

784 1. Although the poleward boundary for the ionospheric potential is specified at the 

785 projection of geostationary orbit in most models (probably except the Rice Convection 

786 Model), we are not able to specify well the ionospheric potential there. Indeed, the 

787 existing models of ionospheric electric potential (like the AMIE technique [ Richmond 

788 and Kamide, 1988], the Weimer [1996, 2001] and the Boyle et al, [1997] models) are 

789 much more reliable at high latitudes and give a poor representation of the potential and 

790 its significant variation in the inner magnetosphere [ Foster and Vo, 2002], In addition, 

791 the effect of FACs is proportional to the volume of the magnetic flux tube, and so 

792 this effect at L=6.6 is about 20% of the FAC effect at L=10, suggesting that FACs 

793 beyond geostationary orbit may produce a major shielding of midlatitudes from a high 

794 latitude driving field. So the region beyond geostationary orbit should be included in 

795 the magnetosphere-ionosphere coupling. An extension of the simulation domain, at least 

796 to A = 72°, is vital for a truly self-consistent modeling of the magnetosphere-ionosphere 

797 coupling. 

798 2. Compared to the case of the Hardy et al. model, the Ridley et al. empirical 

799 relationship between the FAC and conductance produces quite a bit of difference in 
boo the potential distribution and, overall, stronger convection at the subauroral latitudes 
sol (see Figures 4 and 5). This difference strongly affects the cold plasma distribution, 

802 RC precipitation pattern, and EMIC waves (see Figures 7, 11, 12, 13, and 14). More 
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803 importantly, a self-consistent description of the ionospheric conductance makes equation 

804 ( 9 ) nonlinear compared to the case of a statistical conductance model. This is a principle 

805 point requiring that a self-consistent model, based on first principles, of the ionospheric 

806 conductance should be incorporated into a simulation of the magnetosphere-ionosphere 

807 coupling. 

80s 3 . A fine density structure in the plasmasphere boundary layer, plume, detached 

809 plasma etc. controls the wave propagation. This fine structure may be a more crucial 
sio factor in controlling the generation of EMIC waves, than just the intensity/distribution 
an of the RC and the local plasma density. There is very large difference between the wave 

812 activity in the second and third rows in Figures 14 while the density distributions in 

8 1 3 the second and third rows in Figures 7 do not differ so dramatically. This suggests 

aw that to model the EMIC wave distribution and wave spectral properties accurately, the 
sis plasmasphere should be simulated self-consistently because its fine structure requires as 
8 ie much care as that of the RC. 

8 i 7 4 . It is shown that the effect of a finite time needed to reestablish a new potential 

sis pattern throughout the ionosphere and to communicate between the ionosphere and 

8 1 9 the equatorial magnetosphere is important. This effect was ignored in all previous 

820 simulations but it should be taken into account to model a self-consistent electric field 

821 properly. 

822 Concluding we would like to emphasize that in order to make significant progress 

823 in developing a truly self-consistent model of the electric field, we need to considerably 

824 improve our ability to accurately specify the electric field at high latitudes and 
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825 ionospheric conductance. Without this ability, we will not be able to accurately specify 

826 EMIC wave spectra in the inner magnetosphere and correctly describe the wave-induced 

827 heating and/or scattering of the magnetospheric particles. 
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Figure 1. The block diagram of the RC, EMIC waves, plasmasphere, and ionosphere 

H 82 coupling in our model. The system characteristics in orange boxes are externally specified 

and the dashed lines connect the model elements that are currently not linked. 

Figure 2. The interplanetary and geomagnetic characteristics during May 2-4, 1998. 

From the top to the bottom panels: the interplanetary magnetic field GSM By and Bz 

1183 

components, the solar wind velocity, 3-hour Kp index, and the measured Dst index. The 
hours shown are counted from 0000 UT on 1 May, 1998. 

Figure 3. The cross polar cap potential drop from differently driven convection models 
during May 2-4, 1998. The black line, shown for reference, is the potential drop from the 
shielded Volland-Stern model with Kp parameterization. The red, green, and blue lines 
represent the self-consistent results obtained with either the VS or W96 model imposed 

1184 

at A = 69°, and either the Hardy et al. conductance model or the Ridley et al. empirical 
relationship between the FAC and conductance (see legend in the figure). In order to 
drive the W96 model, a 30 min time lag between WIND and the high latitude ionospheric 
boundary is adopted after Farrugia et al. [2003]. 
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Figure 4. The equatorial potential contours in the inner magnetosphere without coro- 
tation held. The view is over the North Pole with local noon to the left. All of the 
indicated hours are counted from 0000 UT on 1 May, 1998. (first row) Results from a 
simulation with the VS model at the high latitude ionospheric boundary and the Hardy et 

1185 

al. conductance model, (second row) Simulation with the W96 model at A = 69° and the 
Hardy et al. conductance model, (third row) The same as in the second row except that 
the Ridley et al. empirical relationship between the FAC and the local Hall/Pedersen 
conductance is used. Equipotentials are drawn every 8 kV. 

Figure 5. (a, b) The potential profiles on the dawn-dusk meridian, and (c, d) the 

1186 

equatorial radial electric field along MLT=18 for hours 33 and 77. 

H87 Figure 6. Same as Figure 4, except that the corotation field is included. 

Figure 7. The equatorial cold plasma density distributions from three self-consistent 

simulations, (first row) Results from a simulation with the VS model at the high latitude 
ionospheric boundary and the Hardy et al. conductance model, (second row) Simulation 

1188 

with the W96 model at A = 69° and the Hardy et al. conductance model, (third row) The 
same as in the second row except that the Ridley et al. empirical relationship between 
the FAC and conductance is used. 

Figure 8. The equatorial cold plasma density distribution in the extended domain of 
L < 10. The electric field is specified by the W96 model above A = 69° but it is calculated 

1189 

self-consistent ly below this latitude using the Ridley et al. relationship between the FAC 


and conductance. 
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Figure 9. The total radial electric field (including the corotation field) in the equatorial 
plane. A combination of the W96 model and the Ridley et al. relationship was used to 

1190 produce these results. Two profiles for MLT=18 and 19 are shown for hours 28 and 29. 
The positive (negative) radial electric field is considered to be parallel (antiparallel) to 
the radius- vector. 

Figure 10. The total equatorial radial electric field versus MLT. A combination of the 
W96 model and the Ridley et al. relationship was used to produce these results. Three 

1191 

profiles for L=8, 9, and 10 are shown for hour 77. The positive (negative) radial electric 

field is considered to be parallel (antiparallel) to the radius-vector. 

Figure 11. The equatorial cold plasma density versus L-shcll for hours 33 and 77. The 

H92 profiles for hour 33 are plotted along MLT=19, while the profiles for hour 77 are plotted 

along MLT=18. 

Figure 12. The RC proton precipitating fluxes averaged over the equatorial pitch-angle 

1193 

loss cone and integrated over the energy range 1 — 50 keV. 

Figure 13. Same as Figure 12, except that the precipitating fluxes are integrated over 

1194 

the energy range 50 — 400 keV. 

Figure 14. The distributions of squared wave magnetic field for the He + - mode EMIC 
waves, (first row) Results from a simulation with the VS model at the high latitude 
ionospheric boundary and the Hardy et al. conductance model, (second row) Simulation 

1195 

with the W96 model at the ionospheric boundary and the Hardy et al. conductance 
model, (third row) The same as in the second row except that the Ridley et al. empirical 
relationship between the FAC and conductance is used. 



60 


1196 


1197 


Figure 15. The equatorial potential contours in the inner magnetosphere without a 
corotation field. The view is over the North Pole with local noon to the left. All of 
the results are from simulations with the W96 potential at the high latitude ionospheric 
boundary and use the Hardy et al. conductance model, (first row) The magnetospheric 
electric field is updated each minute in accordance with the instantaneous interplanetary 
conditions (a 30 min time delay is applied) and FACs, (second row) The interplanetary 
parameters and FACs are averaged over a 20 min window prior to sending them to the 
ionospheric solver and the magnetospheric electric field is updated once every 20 min. 
Equipotentials are drawn every 8 kV. 

Figure 16. The equatorial cold plasma density distributions from simulations with the 
W96 potential at the high latitude ionospheric boundary and the Hardy et al. con- 
ductance model, (first row) The magnetospheric electric field is updated each minute 
accordingly to the instantaneous interplanetary conditions (with a 30 min time delay) 
and FACs, (second row) The interplanetary parameters and FACs are averaged over a 
20 min window prior to sending them to the ionospheric solver and the magnetospheric 


electric field is updated once every 20 min. 
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